02A_MetaPhlanTaxonomic2AMP

library(shades)

Attaching package: ‘shades’

The following object is masked from ‘package:phyloseq’:

    distance

The following object is masked from ‘package:ape’:

    complement

Read in and Prep-process data

Read in Taxa Tables

taxa_table <- read_table('../../results/metaphlan_bowtie_out/merged_metaphlan_profile.tsv',skip = 1) %>%
                          select(!clade_taxid)


# Remove .metaphlan_profile from ends of files 
colnames(taxa_table) <- sapply(strsplit(names(taxa_table), ".metaphlan_profile"), `[[`, 1)
taxa_table <- taxa_table %>% select(!"20221221-Comunal-Zymo")
genus_table <- taxa_table %>% filter(grepl("g__", clade_name), !grepl("s__", clade_name))
genus_table <- genus_table %>% column_to_rownames(var = c("clade_name"))

species_table <- taxa_table %>% filter(grepl("s__", clade_name))
species_table <- species_table %>% column_to_rownames(var = c("clade_name"))

find samples with 100% unknown and remove them

samples.unknown <- species_table %>% 
  melt () %>% 
  group_by(variable) %>% 
  summarise(sum=sum(value)) %>% 
  filter(sum == 0) %>% 
  pluck(1) %>% 
  as.vector()

# Remove samples with no known taxa
species_table <- species_table[, -which(names(species_table) %in% samples.unknown)] 

read in map table

mgx <- read_csv("../../results/minimal_SCRIPT_01_Metagenomics.csv", name_repair = "universal") %>%
  select(Sample, DNA.tube.ID, Batch_sample) %>%
  rename(tube = "DNA.tube.ID", file_name = "Batch_sample") %>%
  filter(file_name %in% colnames(species_table)) %>%
  mutate(Sample = stringr::str_replace_all(Sample, "-", ""))
amp <- clusters %>%
  separate(name, into=c("Sample", "tube"), sep = "-TUBE") %>%
  mutate(Sample = stringr::str_replace_all(Sample, "PT", "")) %>%
  mutate(cluster_num = as.factor(cluster_num)) %>%
  filter(tube %in% mgx$tube)
map_table <- mgx %>% 
  left_join(amp) %>%
  column_to_rownames("file_name") %>% 
  filter(!rownames(.) %in% samples.unknown) %>%
  filter(!is.na(cluster_num))

remove neg controls; set sample.type factor

species_table <- species_table[,rownames(map_table)] # reorder to match metadata

genus_table <- genus_table[,rownames(map_table)]

all(row.names(map_table) %in% colnames(species_table))

Testing as phyloseq + Playground

MPA

library(tidyverse); packageVersion("tidyverse")     #version:1.3.0 
library(phyloseq); packageVersion("phyloseq")       #version:1.32.0


#Building MetaPhlAn species abundance ps object
s_abund <- species_table %>%
  rownames_to_column(var="taxonomy")

s_tax_tab <- s_abund %>%
  dplyr::select(taxonomy) %>%
  tidyr::separate(taxonomy, into = c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species"), sep = "\\|") %>%
  dplyr::mutate(spec_row = Species) %>%
  tibble::column_to_rownames(var = "spec_row") %>%
  filter(Kingdom == "k__Bacteria")

s_otu_tab <- s_abund %>%
  tidyr::separate(taxonomy, into = c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species"), sep = "\\|") %>%
  filter(Kingdom != "k__Bacteria") %>%
  dplyr::select(-Kingdom, -Phylum, -Class, -Order, -Family, -Genus) %>%
  tibble::column_to_rownames(var == "Species")

names(s_otu_tab) <- gsub(names(s_otu_tab), pattern = ".metaphlan_profile", replacement = "") 

head(colSums(s_otu_tab))
s_otu_tab <- s_otu_tab / 100                                                   #convert to proportion with unit sum of 1
head(colSums(s_otu_tab))

s_meta <- data.frame(seq_id = names(s_otu_tab))
s_meta <- s_meta %>%
  left_join(
   rownames_to_column(map_table), 
    by=c("seq_id" = "rowname")) %>%  
  dplyr::mutate(sampleNames_row = seq_id) %>%
  tibble::column_to_rownames(var = "sampleNames_row") %>%
  mutate(sample.id = paste0("PT", Sample, "-TUBE", tube))


(ps_mpa3 <- phyloseq(sample_data(s_meta),
                             otu_table(s_otu_tab, taxa_are_rows = TRUE),
                             tax_table(as.matrix(s_tax_tab))))
ps_mpa3 <- prune_samples(sample_sums(ps_mpa3) > 0, ps_mpa3)
ps_mpa3 <- microbiome::transform(ps_mpa3,  
                                 transform = "compositional",
                                 target="Genus")
ps_mpa3.gen <- tax_glom(ps_mpa3, taxrank = "Genus",NArm=FALSE)

library(vegan)
decostand(s_otu_tab, method="normalize", MARGIN=1)
View(otu_table(ps_mpa3))

KRAKEN

Bray Curtis Distance

dist.amp.bray <-as.matrix(distance(subset_samples(ps.gen.comp, sample.id %in% s_meta$sample.id), method="bray", type="samples"))
dist.mpa.bray <-as.matrix(distance(ps_mpa3.gen, method="bray", type="samples"))
dist.krak.bray <-as.matrix(distance(ps_kraken_species.css, method="bray", type="samples"))

PCoA

PCoA

jacc.pcoa <- pcoa(dist.amp.bray)
pcoa.pt <- jacc.pcoa$vectors %>%
  as_tibble(rownames = NA) %>% 
  rownames_to_column()

pcoa.pt <- pcoa.pt %>% select(1:20) %>%
  left_join(s_meta, by=c("rowname"="sample.id"))

p.wu <-ggplot(pcoa.pt, aes(Axis.1, Axis.2))+ 
    geom_point(aes(color=factor(cluster_num)),size=1, alpha=0) + 
    theme_bw(10) + 
    theme(legend.position = "left") + 
    labs(x="PCoA1",y="PCoA2") +
    #scale_color_brewer(palette="Set1")+  
    scale_color_npg() +
    theme(axis.title = element_text(size = 8)) + 
    geom_point(aes(fill=factor(cluster_num)),size=3, alpha=1, shape=21)  + 
    scale_fill_npg() +
    #scale_fill_brewer(palette="Set1") + 
    nature_theme("PCoA1", "PCoA1") + 
    theme(legend.position = "bottom") + coord_fixed(ratio=1) + ggtitle("Amplicon Bray Curtis")

p.wu

jacc.pcoa <- pcoa(dist.mpa.bray)
pcoa.pt <- jacc.pcoa$vectors %>%
  as_tibble(rownames = NA) %>% 
  rownames_to_column()

pcoa.pt <- pcoa.pt %>% select(1:20) %>%
  left_join(s_meta, by=c("rowname"="seq_id"))

Error in `select()`:
! Can't subset columns past the end.
ℹ Locations 9, 10, 11, …, 19, and 20 don't exist.
ℹ There are only 8 columns.
Backtrace:
 1. pcoa.pt %>% select(1:20) %>% ...
 4. dplyr:::select.data.frame(., 1:20)

UMAP UniFrac W

# Unifrac Settings 
library(ggExtra)
library(umap)

# Set UMAP config (see knight 2021 umap paper)
umap.config<-umap::umap.defaults
umap.config$n_neighbors <- length(sample_names(ps_kraken_species))-1
umap.config$min_dist <- .9
umap.config$spread <- 1
umap.config$random_state <- 1234567
umap.config$transform_state <- 123

do_umap_mgx2amp <- function(dist.matrix, umap_config,join_var){
  # Run UMAP
  fun.umap <- umap::umap(dist.matrix, config = umap_config,input="dist")
  
  # Get Coordinates for ggplot2
  umap.pt <- fun.umap$layout %>%  
    as_tibble(rownames = NA) %>% 
    rownames_to_column()

  # Join with extra metadata, clusters, etc.
  umap.pt <- umap.pt %>% 
    left_join(s_meta, by=c("rowname"=join_var))

  
  # Plot UMAP
  p.umap <-ggplot(umap.pt, aes(x=V1, y=V2))+ 
    geom_point(aes(color=factor(cluster_num)),size=1, alpha=0) + 
    theme_bw(10) + 
    theme(legend.position = "left") + 
    labs(x="UMAP1",y="UMAP2") +
    #scale_color_brewer(palette="Set1")+  
    scale_color_npg() +
    theme(axis.title = element_text(size = 8)) + 
    geom_point(aes(fill=factor(cluster_num)),size=2, alpha=1, shape=21)  + 
    scale_fill_npg() +
    #scale_fill_brewer(palette="Set1") + 
    nature_theme("UMAP1", "UMAP2") + 
    theme(legend.position = "bottom") + theme(aspect.ratio=1)
  p.ellipse <- p.umap + stat_ellipse(aes(color = factor(cluster_num))) #+ theme(legend.position = "none")
  # Optional Plot with Marginal Distributions
  p.marginal<-ggMarginal(p.umap, type="density", groupFill = T, groupColour = T)
  return(list(p.umap, p.ellipse, p.marginal, umap.pt))
}

p.umaps.amp <- do_umap_mgx2amp(dist.amp.bray, umap.config, "sample.id")
p.umaps.amp[[3]]

p.umaps.krak <- do_umap_mgx2amp(dist.krak.bray, umap.config, "seq_id")
p.umaps.krak[[4]] %>%
  left_join(meta) %>%
  ggplot(., aes(x=V1, y=V2))+ 
    geom_point(aes(color=factor(cluster_num)),size=1, alpha=0) + 
    geom_point(aes(fill=factor(cluster_num)),size=2, alpha=1, shape=21)  + 
    stat_ellipse(aes(color = factor(cluster_num))) + 
    scale_color_npg() +
    scale_fill_npg() +
    labs(x="UMAP1",y="UMAP2") +
    nature_theme("UMAP1", "UMAP2") + 
    theme(aspect.ratio=1,legend.position = "bottom", axis.title = element_text(size = 8))
p.umaps.krak[[3]]
p.umaps.krak[[4]] %>%
  left_join(meta) %>%
  ggplot(., aes(x=V1, y=V2))+ 
    #geom_point(color="grey", size=1, alpha=.01) + 
    geom_point(aes(fill=antibiotics_48h),size=2, alpha=1, shape=21)  + 
    scale_fill_npg() +
    labs(x="UMAP1",y="UMAP2") +
    nature_theme("UMAP1", "UMAP2") + 
    theme(aspect.ratio=1,legend.position = "bottom", axis.title = element_text(size = 8))

library(ade4)

tmp.mant <- mantel.rtest(as.dist(dist.amp.bray), as.dist(dist.krak.bray), nrepet = 999)
tmp.mant$^2

Visualize Specific Features

psmelt(ps_kraken_species) %>% 
  arrange(Sample) %>%
  mutate(OTU = stringr::str_replace(string = OTU, pattern = "-", replacement=".")) %>% 
  #filter(Sample %in% c("DNA_B02_27")) %>%
  filter(str_detect(OTU, "k__Viruses")) %>%
  group_by(cluster_num) %>%
  summarize(min = min(Abundance),
            q1 = quantile(value, 0.25),
            median = median(value),
            mean = mean(value),
            q3 = quantile(value, 0.75),
            max = max(value))

Warning in psmelt(ps_kraken_species) :
  The sample variables: 
Sample
 have been renamed to: 
sample_Sample
to avoid conflicts with special phyloseq plot attribute names.

Error in `summarize()`:
ℹ In argument: `q1 = quantile(value, 0.25)`.
ℹ In group 1: `cluster_num = 1`.
Caused by error in `quantile()`:
! object 'value' not found
Backtrace:
  1. ... %>% ...
 11. stats::quantile(value, 0.25)

Maaslin Genera

meta_4mas <- sample_data(ps_kraken_species) %>%
  as.data.frame() %>%
  as_tibble(rownames = NA) %>%
  select(sample.id, seq_id, cluster_num) %>% 
  left_join(select(meta, sample.id, pt_category, success,patient.id.num)) %>%
  column_to_rownames("seq_id") #%>% mutate(cluster_num = if_else(cluster_num == "4", "0", cluster_num))

Joining with `by = join_by(sample.id)`

krakenMGX = Maaslin2(input_data     = kraken_4mas, 
                     input_metadata = meta_4mas, 
                     min_abundance  = .001,
                     min_prevalence = .1,
                     #normalization  = "TSS",
                     transform = "AST",
                     output         = "../../results/notebook_out/maaslin2_out", 
                     fixed_effects  = c("cluster_num"),
                     reference      = c("cluster_num,4"),
                     random_effects = c("patient.id.num"),
                     core           = 4,plot_scatter = FALSE
                     ) # Subject

[1] "Warning: Deleting existing log file: ../../results/notebook_out/maaslin2_out/maaslin2.log"
2023-04-26 15:00:14 INFO::Writing function arguments to log file
2023-04-26 15:00:14 INFO::Verifying options selected are valid
2023-04-26 15:00:14 INFO::Determining format of input files
2023-04-26 15:00:14 INFO::Input format is data samples as rows and metadata samples as rows
2023-04-26 15:00:14 INFO::Formula for random effects: expr ~ (1 | patient.id.num)
2023-04-26 15:00:14 INFO::Formula for fixed effects: expr ~  cluster_num
2023-04-26 15:00:14 INFO::Filter data based on min abundance and min prevalence
2023-04-26 15:00:14 INFO::Total samples in data: 202
2023-04-26 15:00:14 INFO::Min samples required with min abundance for a feature not to be filtered: 20.200000
2023-04-26 15:00:14 INFO::Total filtered features: 1402
2023-04-26 15:00:14 INFO::Filtered feature names from abundance and prevalence filtering:
2023-04-26 15:00:14 INFO::Total filtered features with variance filtering: 0
2023-04-26 15:00:14 INFO::Filtered feature names from variance filtering:
2023-04-26 15:00:14 INFO::Running selected normalization method: TSS
2023-04-26 15:00:14 INFO::Applying z-score to standardize continuous metadata
2023-04-26 15:00:14 INFO::Running selected transform method: AST
2023-04-26 15:00:14 INFO::Running selected analysis method: LM
2023-04-26 15:00:14 INFO::Creating cluster of 4 R processes

  |                                                  | 0 % ~calculating  
  |+++                                               | 6 % ~43s          
  |++++++                                            | 11% ~21s          
  |+++++++++                                         | 17% ~14s          
  |++++++++++++                                      | 22% ~11s          
  |++++++++++++++                                    | 28% ~08s          
  |+++++++++++++++++                                 | 33% ~07s          
  |++++++++++++++++++++                              | 39% ~06s          
  |+++++++++++++++++++++++                           | 44% ~05s          
  |+++++++++++++++++++++++++                         | 50% ~04s          
  |++++++++++++++++++++++++++++                      | 56% ~03s          
  |+++++++++++++++++++++++++++++++                   | 61% ~03s          
  |++++++++++++++++++++++++++++++++++                | 67% ~02s          
  |+++++++++++++++++++++++++++++++++++++             | 72% ~02s          
  |+++++++++++++++++++++++++++++++++++++++           | 78% ~01s          
  |++++++++++++++++++++++++++++++++++++++++++        | 83% ~01s          
  |+++++++++++++++++++++++++++++++++++++++++++++     | 89% ~01s          
  |++++++++++++++++++++++++++++++++++++++++++++++++  | 94% ~00s          
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=06s  
2023-04-26 15:00:21 INFO::Counting total values for each feature
2023-04-26 15:00:21 WARNING::Deleting existing residuals file: ../../results/notebook_out/maaslin2_out/residuals.rds
2023-04-26 15:00:21 INFO::Writing residuals to file ../../results/notebook_out/maaslin2_out/residuals.rds
2023-04-26 15:00:21 WARNING::Deleting existing fitted file: ../../results/notebook_out/maaslin2_out/fitted.rds
2023-04-26 15:00:21 INFO::Writing fitted values to file ../../results/notebook_out/maaslin2_out/fitted.rds
2023-04-26 15:00:21 WARNING::Deleting existing ranef file: ../../results/notebook_out/maaslin2_out/ranef.rds
2023-04-26 15:00:21 INFO::Writing extracted random effects to file ../../results/notebook_out/maaslin2_out/ranef.rds
2023-04-26 15:00:21 INFO::Writing all results to file (ordered by increasing q-values): ../../results/notebook_out/maaslin2_out/all_results.tsv
2023-04-26 15:00:21 INFO::Writing the significant results (those which are less than or equal to the threshold of 0.250000 ) to file (ordered by increasing q-values): ../../results/notebook_out/maaslin2_out/significant_results.tsv
2023-04-26 15:00:21 INFO::Writing heatmap of significant results to file: ../../results/notebook_out/maaslin2_out/heatmap.pdf

Warning in xtfrm.data.frame(x) : cannot xtfrm data frames

maaslin_norm.gen <- krakenMGX$results %>% as.data.frame() 
maaslin_norm.gen$OTU <- rownames(maaslin_norm.gen)
maaslin_norm.gen <- maaslin_norm.gen %>% 
  #separate(feature, into=c("ASV", "feature"), sep = "[.]") %>% 
  filter(qval < 0.00001) %>%
  filter(N.not.zero > 10) %>%
  arrange(qval) %>% 
  select(feature) %>% 
  unique
diff.genera <- maaslin_norm.gen$feature
sort(diff.genera)

 [1] "Corynebacterium_genitalium"       "Corynebacterium_kefirresidentii" 
 [3] "Corynebacterium_singulare"        "Corynebacterium_striatum"        
 [5] "Cutibacterium_acnes"              "Gemella_haemolysans"             
 [7] "Lacticaseibacillus_rhamnosus"     "Lancefieldella_parvula"          
 [9] "Mogibacterium_diversum"           "Moraxella_osloensis"             
[11] "Olsenella_uli"                    "Parvimonas_micra"                
[13] "Prevotella_jejuni"                "Rothia_mucilaginosa"             
[15] "Schaalia_odontolytica"            "Staphylococcus_capitis"          
[17] "Staphylococcus_epidermidis"       "Staphylococcus_hominis"          
[19] "Streptococcus_anginosus"          "Streptococcus_parasanguinis"     
[21] "Streptococcus_salivarius"         "Streptococcus_sp._HSISM1"        
[23] "Streptococcus_sp._LPB0220"        "Streptococcus_sp._oral_taxon_061"
[25] "Streptococcus_sp._oral_taxon_431"

otu.ast.krak.phy <- psmelt(ps.ast.krak.phy) %>% 
  mutate(Abundance = sign(Abundance)*asin(sqrt(abs(Abundance)))) %>% # if want AST transform
  arrange(Sample) %>%
  select(Sample, OTU, Abundance) %>%
  pivot_wider(id_cols = OTU, names_from = Sample, values_from = Abundance) %>%
  column_to_rownames(var="OTU") %>% 
  as.matrix() 

Warning in psmelt(ps.ast.krak.phy) :
  The sample variables: 
Sample
 have been renamed to: 
sample_Sample
to avoid conflicts with special phyloseq plot attribute names.

psmelt(ps_kraken_phyla) %>% 
  select(Sample,Abundance, Species,cluster_num) %>%
# mutate(k__Viruses = log10(k__Viruses),
  #        k__Bacteria = log10(k__Bacteria),
  #        k__Eukaryota = log10(k__Eukaryota)) %>%
  ggplot(., aes(x=cluster_num, y=Abundance, fill=Species, group=cluster_num)) +
  geom_boxplot(aes(fill=Species), outlier.alpha = 0) +
  geom_point(aes(fill=Species), 
             alpha = 0.75,
             size = 1, 
             shape = 21, 
             color = 'black',
             position = position_jitterdodge(jitter.width = .5)) +
  scale_fill_simpsons() +
  nature_theme("Read", "Fol") +
  labs(x="Pneumotype Cluster" , 
       y="Domain Abundance", 
       title="Domain Level Shifts in Metagenomic Data  (AST Normalized)") +
  theme(aspect.ratio = 1) + facet_grid(~Species) + 
  stat_compare_means(label.y.npc = .95,size=2,label.x.npc = "center")

Warning in psmelt(ps_kraken_phyla) :
  The sample variables: 
Sample
 have been renamed to: 
sample_Sample
to avoid conflicts with special phyloseq plot attribute names.